ACCESS overestimates concentration in polar region while underestimates in
import cosima_cookbook as cc
from cosima_cookbook import explore
import numpy as np
import dask.array as da
import xarray as xr
import pandas as pd
import os.path
import glob
from tqdm.notebook import tqdm
#------------------------
#----------------------
import logging
logging.captureWarnings(True)
logging.getLogger('py.warnings').setLevel(logging.ERROR)
from dask.distributed import Client
# change to your own directory on /g/data/:
figdir = '/g/data/hh5/tmp/rsp599_tmp/ACESSOM201/figures/BGC_IAF/'
Recommended plotting libraries
%matplotlib inline
%config InlineBackend.figure_format='retina'
#Packages for plotting
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.dates as mdates
import matplotlib.colors as mcolors
from matplotlib.ticker import AutoMinorLocator
import matplotlib.cm as mcm
import matplotlib.gridspec as gridspec
#
import IPython.display
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import geopandas
#
land_50m = cft.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='black', facecolor='gray', linewidth=0.5)
client = Client()
client
# client = Client(processes=2, threads_per_worker=1,
# n_workers=2, memory_limit='10GB')
# client
Client-7341e282-86a0-11ed-9493-fa163eefa180
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: /proxy/8787/status |
be1fde5a
| Dashboard: /proxy/8787/status | Workers: 4 |
| Total threads: 16 | Total memory: 44.92 GiB |
| Status: running | Using processes: True |
Scheduler-06c0426a-5aff-474e-a57b-f603dd683825
| Comm: tcp://127.0.0.1:36025 | Workers: 4 |
| Dashboard: /proxy/8787/status | Total threads: 16 |
| Started: Just now | Total memory: 44.92 GiB |
| Comm: tcp://127.0.0.1:41483 | Total threads: 4 |
| Dashboard: /proxy/46387/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:46659 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-wtxqdkwu | |
| Comm: tcp://127.0.0.1:37431 | Total threads: 4 |
| Dashboard: /proxy/37311/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:46353 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-pv5unken | |
| Comm: tcp://127.0.0.1:36463 | Total threads: 4 |
| Dashboard: /proxy/34123/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:36187 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-ve5kdln6 | |
| Comm: tcp://127.0.0.1:38027 | Total threads: 4 |
| Dashboard: /proxy/42683/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:38931 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-2jcanrz5 | |
Wilma has kindly prepared the world ocean atlas data, currently 0.1$^\circ$ version available at /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/
!ncdump -h /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o01_10.nc
netcdf woa18_all_o01_10 {
dimensions:
xt_ocean = 3600 ;
yt_ocean = 2700 ;
st_ocean = 51 ;
time = 1 ;
variables:
double xt_ocean(xt_ocean) ;
xt_ocean:_FillValue = NaN ;
xt_ocean:long_name = "Nominal Longitude of T-cell center" ;
xt_ocean:units = "degree_east" ;
xt_ocean:modulo = 360. ;
xt_ocean:point_spacing = "even" ;
xt_ocean:axis = "X" ;
double yt_ocean(yt_ocean) ;
yt_ocean:_FillValue = NaN ;
yt_ocean:long_name = "tcell latitude" ;
yt_ocean:units = "degrees_N" ;
yt_ocean:cartesian_axis = "Y" ;
double st_ocean(st_ocean) ;
st_ocean:_FillValue = NaN ;
st_ocean:long_name = "zt" ;
st_ocean:units = "meters" ;
st_ocean:positive = "down" ;
st_ocean:point_spacing = "uneven" ;
st_ocean:axis = "Z" ;
double time(time) ;
time:_FillValue = NaN ;
time:long_name = "time" ;
time:cartesian_axis = "T" ;
time:calendar_type = "GREGORIAN" ;
time:units = "days since 0001-01-01" ;
time:calendar = "GREGORIAN" ;
double o_an(time, st_ocean, yt_ocean, xt_ocean) ;
o_an:_FillValue = NaN ;
o_an:long_name = "Objectively analyzed mean fields for mole_concentration_of_dissolved_molecular_oxygen_in_sea_water at standard depth levels." ;
o_an:units = "micromoles_per_kilogram" ;
}
oxy = xr.open_mfdataset('/g/data/hh5/tmp/rsp599_tmp/oxygen_srf_mon_clim_2000_2019.nc')
oxy
<xarray.Dataset>
Dimensions: (xt_ocean: 3600, yt_ocean: 2700, month: 12)
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
* month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
Data variables:
Oxygen (month, yt_ocean, xt_ocean) float32 dask.array<chunksize=(12, 2700, 3600), meta=np.ndarray>
Attributes:
long_name: dissolved oxygen
units: mmol/m^3
valid_range: [-1000000. 1000000.]
cell_methods: time: mean
time_avg_info: average_T1,average_T2...
QuantizeGranularBitRoundNumberOfSignificantDigits: 4
ncfiles: ['/g/data/cj50/access...
contact: Andrew Kiss
email: andrew.kiss@anu.edu.au
created: 2022-04-27
description: 0.1 degree ACCESS-OM2...
notes: Run configuration and...woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o01_10.nc')
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
<xarray.DataArray 'o_an' (st_ocean: 51, yt_ocean: 2700, xt_ocean: 3600)>
dask.array<getitem, shape=(51, 2700, 3600), dtype=float64, chunksize=(51, 2700, 3600), chunktype=numpy.ndarray>
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
* st_ocean (st_ocean) float64 0.5413 1.681 2.94 ... 1.333e+03 1.453e+03
Attributes:
long_name: Objectively analyzed mean fields for mole_concentration_of_di...
units: micromoles_per_kilogramoxydiff = oxy.Oxygen.isel(month=0) - woa_oxy.o_an.isel(st_ocean=0)
oxydiff.load()
<xarray.DataArray (yt_ocean: 2700, xt_ocean: 3600, time: 1)>
array([[[ nan],
[ nan],
[ nan],
...,
[ nan],
[ nan],
[ nan]],
[[ nan],
[ nan],
[ nan],
...,
[ nan],
[ nan],
[ nan]],
[[ nan],
[ nan],
[ nan],
...,
...
...,
[25.60482412],
[25.44903289],
[25.2837304 ]],
[[25.39055891],
[25.22703201],
[25.05789347],
...,
[25.88977886],
[25.72440306],
[25.5567456 ]],
[[25.67717358],
[25.51730854],
[25.35306159],
...,
[26.17662515],
[26.01076811],
[25.84190407]]])
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
month int64 1
st_ocean float64 0.5413
* time (time) object 0001-01-02 00:00:00fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Jan Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o02_10.nc')
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=1) - woa_oxy.o_an.isel(st_ocean=0)
oxydiff.load()
#
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Feb Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o03_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=2) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('March Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o04_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=3) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('April Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o05_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=4) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('May Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o06_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=5) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('June Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o07_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=6) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('July Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o08_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=7) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Aug Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o09_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=8) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Sep Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o10_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=9) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Oct Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o11_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=10) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Nov Clim (ACCES - WOA)');
woa_oxy = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o12_10.nc') # change filename
woa_oxy = woa_oxy.assign_coords(xt_ocean=(((woa_oxy.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_oxy.o_an.squeeze(drop=True)
#
oxydiff = oxy.Oxygen.isel(month=11) - woa_oxy.o_an.isel(st_ocean=0)### Change month name
oxydiff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$Oxygen (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
oxydiff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.balance, cbar_kwargs=cbar_kwargs)
plt.title('Dec Clim (ACCES - WOA)');